/**
  SUNNY Project, Anyang Normal University, IMP-CAS
  \class TAHamiltonian
  \brief The Hamiltonian matrix class, to generate and manage the many-body Slater
  determinants, and also interfaces the matrix element calculation methods. Note
  that this is an abstract base class. Calculation of the matrix elements should
  be specifically implemented in its sub-class.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/02/09
  \date Last modified: 2021/02/04, by SUN Yazhou
  \copyright 2020-2021 SUN Yazhou
  \copyright SUNNY project, Anyang Normal University, IMP-CAS
*/

#include <algorithm>
#include "TAMBSD.h"
#include "TAHamiltonian.h"
#include "TAMBSDList.h"
#include "TAMBSDManager.h"
#include "TASPState.h"
#include "TAParaManager.h"

using std::swap;

static const double DUMB = -9999.;

TAHamiltonian::TAHamiltonian(const string &configFile)
    : fMBSDListM(0), fNSPState(0), fMatrix(nullptr){
  // prepare the single-particle basis set //
  // prepare the basis of the representation //
  TAMBSDManager *mbsdManager = TAMBSDManager::Instance();
  mbsdManager->LoadConfigFile(PARA->BaseDir()+configFile);
  mbsdManager->MSchemeGo(); // generate many-body basis
  fMBSDListM = mbsdManager->GetMBSDListM();
  fNSPState = mbsdManager->GetNSPState();
  fNParticle = mbsdManager->GetNParticle();
} // end the constructor

TAHamiltonian::~TAHamiltonian(){} // end the destructor

ull TAHamiltonian::GetNBasis() const{
  return fMBSDListM->GetNBasis();
} // end member function GetNBasis

/// calculate H(rr,cc) \param r: row, c: column
double TAHamiltonian::MatrixElement(int rr, int cc) const{
  if(rr < cc) swap(rr, cc);

  return MatrixElement1N(rr, cc) + MatrixElement2N(rr, cc) +
    MatrixElement3N(rr, cc);
} // end member function MatrixElement

// r = H*v
void TAHamiltonian::DotProduct(const TASparseVec &v, TASparseVec &r) const{
  const int n = GetNBasis();
  r.Clear(); // clear all the entries
  double t, e1, e2;
  for(int i = 0; i < n; i++){
    t = 0.;
    for(int j = 0; j < n; j++)
      if((e1 = MatrixElement(i,j)) && (e2 = v[j])) t += e1 * e2;
    r.Fill(i,t);
  } // end for over i
} // end member function DotProduct

/// construct the Hamiltonian matrix
matrix *TAHamiltonian::Matrix(){
  if(fMatrix) return fMatrix;

  const int n = GetNBasis();
  if(n > 10000) TAException::Warn("TAHamiltonian",
    "Matrix: dimension too large (>10000).");
  fMatrix = new matrix(n, n);
  for(int i = 0; i < n; i++) for(int j = 0; j <= i; j++)
    (*fMatrix)[j][i] = (*fMatrix)[i][j] = MatrixElement(i, j);

  return fMatrix;
} // end member function Matrix
